Exploring Age Distributions

Load libraries

library(dplyr)
library(readr)
library(rgdal)
library(spdep)
library(tidyr)
library(ggplot2)
library(geofacet)
library(viridis)
library(plotly)

Read in required files

bairro_shp <- readOGR("Shapefiles/Bairro/", "Bairros_from_CTs")

OGR data source with driver: ESRI Shapefile 
Source: "/Users/Sudipta/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit/Shapefiles/Bairro", layer: "Bairros_from_CTs"
with 119 features
It has 4 fields

Get median age at various levels of stratification (bairro, year, and bairro-year)

#Get overall median age
median_age_all <- median(homs$Age, na.rm = TRUE)
#Create dataframes of median by different strata
age_by_yr_bairro <- homs %>% dplyr::select(Age, Bairro, YOD) %>% group_by(Bairro, YOD) %>% summarize(median_age = median(Age, na.rm=TRUE), count = n(), mean_age = mean(Age, na.rm=TRUE))
age_by_bairro <- homs %>% dplyr::select(Age, Bairro, YOD) %>% group_by(Bairro) %>% summarize(median_age = median(Age, na.rm=TRUE), count = n(), mean_age = mean(Age, na.rm=TRUE))
age_by_yr <- homs %>% dplyr::select(Age, Bairro, YOD) %>% group_by(YOD) %>% summarize(median_age = median(Age, na.rm=TRUE), count = n(), mean_age = mean(Age, na.rm=TRUE))

Now plotting yearly median age and count on the same plot

ggplot(age_by_yr, aes(x=YOD)) + geom_line(aes(y=median_age, colour="Median Age")) + geom_line(aes(y=log(count)*3.9, colour="Count")) + #transform to get on same axis
scale_y_continuous(sec.axis = sec_axis(~exp(./3.9), name="Count")) +  scale_colour_manual(values = c("blue", "red")) + labs(y = "Median Age", x = "Year")

Plot the median age by bairro on a map

bairro_shp@data <- cbind(bairro_shp@data, trueCentroids@coords)
bairro_ggplot <- shape_to_ggplot(bairro_shp)

Regions defined for each Polygons
column name ‘id’ is duplicated in the result

ggplot() + geom_polygon(data = bairro_ggplot, aes(x = long, y = lat, group = group, fill=median_age), size=0.1, color="black",) +        
    scale_fill_viridis(trans="reverse") + geom_point(data = bairro_ggplot, aes(x=x,y=y, color=log(count)*4), size=1.5) + scale_color_viridis(option = "magma") +
    theme(line = element_blank(),
          plot.title = element_text(hjust = 0.5),
          axis.text=element_blank(),   
          axis.title=element_blank(),
          legend.text=element_text(size=8),
          legend.title=element_text(size=8),
          panel.background = element_blank()) + coord_equal()

ggplot() + geom_polygon(data = bairro_ggplot, aes(x = long, y = lat, group = group, fill=log(count)*4), size=0.1, color="black") +         
    scale_fill_viridis() +
    theme(line = element_blank(),
          plot.title = element_text(hjust = 0.5),
          axis.text=element_blank(),   
          axis.title=element_blank(),
          legend.text=element_text(size=8),
          legend.title=element_text(size=8),
          panel.background = element_blank()) + coord_equal()

Get deviation for overall median age for all bairros by year

age_by_yr$dev <- age_by_yr$median_age - median_age_all
age_by_yr_bairro$dev <- age_by_yr_bairro$median_age - median_age_all
age_by_bairro$dev <- age_by_bairro$median_age - median_age_all

mygrid <- data.frame(
  name = c("JACARECANGA", "PIRAMBU", "CRISTO REDENTOR", "BARRA DO CEARA", "CAIS DO PORTO", "MEIRELES", "PRAIA DE IRACEMA", "MOURA BRASIL", "VILA VELHA", "JARDIM IRACEMA", "FLORESTA", "ALVARO WEYNE", "CARLITO PAMPLONA", "VARJOTA", "VICENTE PINZON", "MUCURIPE", "ALDEOTA", "CENTRO", "FARIAS BRITO", "QUINTINO CUNHA", "MONTE CASTELO", "VILA ELLERY", "JARDIM GUANABARA", "PRESIDENTE KENNEDY", "OLAVO OLIVEIRA", "SAO GERARDO/ALAGADICO", "DE LOURDES", "PRAIA DO FUTURO I", "GUARARAPES", "COCO", "ANTONIO BEZERRA", "PADRE ANDRADE", "PARQUE ARAXA", "BENFICA", "JOSE BONIFACIO", "JOAQUIM TAVORA", "DIONISIO TORRES", "RODOLFO TEOFILO", "GENIBAU", "JARDIM AMERICA", "LUCIANO CAVALCANTE", "MANUEL DIAS BRANCO", "PRAIA DO FUTURO II", "CIDADE 2000", "SALINAS", "SAO JOAO DO TAUAPE", "FATIMA", "AMADEU FURTADO", "PARQUELANDIA", "PICI", "PAPICU", "CONJUNTO CEARA I", "JARDIM DAS OLIVEIRAS", "AEROLANDIA", "PARQUE MANIBURA", "SAPIRANGA COITE", "EDSON QUEIROZ", "SABIAGUABA", "DOM LUSTOSA", "AUTRAN NUNES", "BELA VISTA", "CIDADE DOS FUNCIONARIOS", "ALTO DA BALANCA", "PARREAO", "BOM FUTURO", "DAMAS", "CONJUNTO CEARA II", "GRANJA LISBOA", "JOAO XXIII", "DIAS MACEDO", "CAJAZEIRAS", "VILA UNIAO", "COUTO FERNANDES", "PAN AMERICANO", "CAMBEBA", "JOSE DE ALENCAR", "CURIO", "LAGOA REDONDA", "HENRIQUE JORGE", "MONTESE", "AEROPORTO", "GRANJA PORTUGAL", "MONDUBIM", "CANINDEZINHO", "PARQUE DOIS IRMAOS", "PASSARE", "PARQUE IRACEMA", "BARROSO", "GUAJERU", "MESSEJANA", "COACU", "PARQUE SAO JOSE", "BOA VISTA", "ITAPERI", "BOM SUCESSO", "SERRINHA", "BOM JARDIM", "SIQUEIRA", "DEMOCRITO ROCHA", "VILA PERY", "JOCKEY CLUBE", "CONJUNTO ESPERANCA", "JANGURUSSU", "ANCURI", "SAO BENTO", "DENDE", "ITAOCA", "JARDIM CEARENSE", "PARANGABA", "PLANALTO AYRTON SENNA", "PARQUE PRESIDENTE VARGAS", "PAUPINA", "PARQUE SANTA ROSA", "MARAPONGA", "PALMEIRAS", "PARQUE SANTA MARIA", "MANOEL SATIRO", "PREFEITO JOSE WALTER", "PEDRAS"), code=c(1:119),
  row = c(1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 11, 11),
  col = c(10, 9, 8, 7, 16, 13, 12, 11, 6, 7, 8, 9, 10, 15, 16, 14, 13, 12, 11, 6, 10, 9, 7, 8, 7, 8, 16, 17, 14, 15, 5, 6, 9, 10, 11, 12, 13, 9, 5, 10, 14, 17, 18, 15, 13, 12, 11, 8, 7, 6, 16, 4, 14, 12, 15, 16, 17, 18, 6, 5, 7, 13, 11, 10, 9, 8, 4, 3, 6, 12, 13, 10, 8, 7, 14, 16, 15, 17, 5, 9, 11, 4, 3, 2, 8, 9, 12, 13, 15, 14, 16, 6, 10, 7, 5, 11, 4, 3, 6, 9, 5, 11, 13, 14, 15, 8, 10, 12, 7, 4, 7, 13, 8, 11, 12, 9, 10, 11, 12),
  stringsAsFactors = FALSE
)
geofacet::grid_preview(mygrid)

Note: You provided a user-specified grid. If this is a
  generally-useful grid, please consider submitting it to become a
  part of the geofacet package. You can do this easily by calling:
  grid_submit(__grid_df_name__)

age_by_yr_bairro$colour <- ifelse(age_by_yr_bairro$dev < 0, "negative","positive")
age_by_yr_bairro <- age_by_yr_bairro[complete.cases(age_by_yr_bairro),]
ggplot(age_by_yr_bairro) + geom_col(aes(YOD, dev, fill=colour)) + 
  scale_fill_manual(values=c(negative="firebrick1",positive="steelblue")) + facet_geo(~Bairro,grid=mygrid, label="name") + theme_void() + theme(strip.background = element_blank(), strip.text.x = element_text(size=4))

Note: You provided a user-specified grid. If this is a generally-useful grid, please consider submitting it to become a
  part of the geofacet package. You can do this easily by calling: grid_submit(__grid_df_name__)

age_by_bairro$colour <- ifelse(age_by_bairro$dev < 0, "negative","positive")
age_by_bairro <- age_by_bairro[complete.cases(age_by_bairro),]
ggplot(age_by_bairro) + geom_col(aes(x=0,y=dev, fill=colour)) +
  scale_fill_manual(values=c(negative="firebrick1",positive="steelblue")) +
  facet_geo(~ Bairro, grid = mygrid, label="name") + theme_void() + theme(strip.background = element_blank(), strip.text.x = element_text(size=2))

Note: You provided a user-specified grid. If this is a generally-useful grid, please consider submitting it to become a
  part of the geofacet package. You can do this easily by calling: grid_submit(__grid_df_name__)

Some scatter plots to explore as well

ggplot(age_by_bairro, aes(x=log(count), y=median_age)) + geom_point() + geom_smooth(method="loess")

p <- ggplot(age_by_yr_bairro, aes(x=log(count), y=median_age)) + geom_point() + geom_smooth(method = "loess")
p

ggplotly(p)

We recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`